 #!/usr/bin/env python
# -*- coding: utf-8 -*-
# @Time    : 2018/1/31 16:01
# @Author  : Aries
# @Site    : 
# @File    : generateLUT.py
# @Software: PyCharm Community Edition
import numpy as np
from collections import namedtuple
from Py6S import *  # need for py6s atmos corr
from config import *
import sensor_OLI

def __build_DEM_lut(aeroProfile, atmosProfile, groundRefl, useBRDF, acqtime, lon, lat):
    """
    : rtype: lut
    """
    LUTPicker = namedtuple('LutPicker', ['toa', 'dem', 'sixsCoeffs'])
    lut = list()
    toa_min = 0
    dem_min = 2
    num_toa_steps = 200
    num_dem_steps = 60
    toa_val = toa_min
    for i in range(num_toa_steps):
        dem_val = dem_min
        for j in range(num_dem_steps):
            # print(dem_val,toa_val)
            lut.append(LUTPicker(toa=toa_val, dem=dem_val, sixsCoeffs=
        sensor_OLI.calc6SAtmosCorr(aeroProfile, atmosProfile, groundRefl, dem_val, toa_val, useBRDF, acqtime, lon,
                                   lat)))
            dem_val = dem_val + 0.1
        toa_val = toa_val + 0.01
    return lut

def _getAtmosCorrLUT(lon, lat, acquisitionTime, useBRDF):
    """
    step1
    :rtype: lut
    :call func: build lut
    """
    lut = list()

    # 设置气溶胶模型
    aeroProfile = None
    if lat > 38.8 and lon > 75.8:
        aeroProfile = AeroProfile.PredefinedType(AeroProfile.Desert)  # 沙漠型气溶胶
    elif lat > 34.8 and lat <= 38.8:
        aeroProfile = AeroProfile.PredefinedType(AeroProfile.Stratospheric)  # 高寒型气溶胶
    # elif near-by-city:
    #	s.aero_profile = AeroProfile.PredefinedType(AeroProfile.Urban) #城市型气溶胶
    else:
        aeroProfile = AeroProfile.PredefinedType(AeroProfile.Continental)  # 默认大陆型气溶胶

    # 设置大气模型
    date = acquisitionTime.strftime('%Y-%m-%d')
    atmosProfile = AtmosProfile.FromLatitudeAndDate(lat, date)  # 利用日期和纬度

    # 设置地面反射
    groundRefl = GroundReflectance.HomogeneousLambertian(Spectra.import_from_usgs(
        "http://speclab.cr.usgs.gov/spectral.lib06/ds231/ASCII/V/russianolive.dw92-4.30728.asc"))  # 利用usgs波谱库

    lut = __build_DEM_lut(aeroProfile, atmosProfile, groundRefl, useBRDF, acquisitionTime, lon, lat)
    return lut

def buildLUT(lonCentre, latCentre, acquisitionTime):
    """
	build lut
	:return:
	"""
    lut = _getAtmosCorrLUT(lonCentre, latCentre, acquisitionTime, True)
    lutarr = np.empty([len(lut), 32])
    for i, lutkeeper in enumerate(lut):
        print(i)
        lutarr[i, 0] = lutkeeper.toa
        lutarr[i, 1] = lutkeeper.dem
        coeffs = lutkeeper.sixsCoeffs.ravel()
        for j in range(len(coeffs)):
            lutarr[i, j+2] = coeffs[j]
    np.savetxt('LC08_L1TP_150033_LUT_BRDF.out', lutarr)
    pass


if __name__ == '__main__':
    header = LC08_20160227 + "/LC08_L1TP_150033_20160227_20170329_01_T1_MTL.txt"
    lonCentre, latCentre, acquisitionTime = sensor_OLI._extractHeaderParameters(header)
    print(lonCentre, latCentre, acquisitionTime)
    buildLUT(lonCentre, latCentre, acquisitionTime)